% plot the g/f curve for an izhikevich neuron.
function gfcurve

b=[];
for neuron = 1:4

    for i = 1:100,a(i)=izhikevich(i,0,neuron,true),end,close all,
    b=[b a'];
end

c=[];
for neuron = 1:4

    for i = 1:100,a(i)=izhikevich(-10*i,0,neuron,false),end,close all,
    c=[c a'];
end
figure (1);
plot(b,'-')
xlabel('g (nS)','FontSize',12);
ylabel ('Firing rate (Hz)','FontSize',12);
title ('Conductance/firing rate curve for neuron','FontSize',14);
legend('Pyr old','Pyr 2008','FS 2008','TC 2008')

figure(2);
plot(10*(1:100),c,'-')
xlabel('I (pA)','FontSize',12);
ylabel ('Firing rate (Hz)','FontSize',12);
title ('Conductance/firing rate curve for neuron','FontSize',14);
legend('Pyr old','Pyr 2008','FS 2008','TC 2008')


% function [spikes Gampa]=izhikevich(stim_conductance_exc, stim_conductance_inh,neuron,g_true_I_false)
% 
% if ~exist('g_true_I_false')
%     g_true_I_false = true;
% end
% if ~exist('neuron')
%     neuron = 1;
% end
% 
% 
% switch neuron
%     case 1
%         C=100;vr=-60;vt=-40;k=3;
%         a=0.01;b=5;c=-50;d=100;
%         vp=35;
%     case 2
%         % pyramidal from Izh & Edelman, 2008
%         vr=-60;
%         vt=-50
%         k=3;
%         C=80
%         a=0.01;
%         b=5;
%         c=-60;
%         d=400;
%         vp=50;
%     case 3
%         % fs cell.
%         vr=-55;
%         vt=-40;
%         k=1;
%         C=20;
%         a=.15;
%         b=8;
%         c=-55;
%         d=200;
%         vp=25;
%     case 4
%         % TC cell from ibid., assuming only tonic firing model
%         vr= -60;
%         vt= -50;
%         k=   1.6;
%         C= 200;
%         a=0.1;
%         b=   0; % Assuming it is in tonic mode all the time.
%         c=   -60;
%         d=   10;
%         vp=  40;
% end
% 
% 
% 
% 
% T=1000;
% 
% dt_ms=1;
% 
% 
% spikes = 0;
% n=round(T/dt_ms);
% v=vr*ones(1,n); u=0*v;
% si= stim_conductance_inh .*[ones(1,n)];
% se= stim_conductance_exc .*[ones(1,n)];
% 
% % model conductances.
% tau_ampa  = exp(-dt_ms/5);
% tau_nmda  = exp(-dt_ms/150);
% tau_gabaa = exp(-dt_ms/6);
% tau_gabab = exp(-dt_ms/150);
% 
% E_AMPA	= 0.0;
% E_NMDA	= 0.0;
% E_GABAa	= -70;
% E_GABAb = -90;
% Ggabaa = 0;Ggabab = 0; Gampa = 0;Gnmda = 0;
% 
% Ggabaa=zeros(1,n);
% Ggabab=zeros(1,n);
% Gampa=zeros(1,n);
% Gnmda=zeros(1,n);
% I=zeros(1,n);
% 
% gain_ampa=1;
% gain_nmda=1;
% gain_gabaa=1;
% gain_gabab=1;
% % end model conductances.
% 
% 
% for i = 1:n-1
%     
%     % conductances from synapses.
%     %Ggabaa(i+1) = Ggabaa(i)*tau_gabaa + dt_ms*(1-tau_gabaa)*gain_gabaa*si(i);
%     %Ggabab(i+1) = Ggabab(i)*tau_gabab + dt_ms*(1-tau_gabab)*gain_gabab*si(i);
%     
%     %Gampa(i+1) = Gampa(i)*tau_ampa + dt_ms*(1-tau_ampa)*gain_ampa*se(i);
%     %Gnmda(i+1) = Gnmda(i)*tau_nmda + dt_ms*(1-tau_nmda)*gain_nmda*se(i);
% 
%     Ggabaa(i+1) = si(i);
%     
%     Gampa(i+1) = se(i);
%     
%     
%     xx = (-80-v(i))/60;
% 	xx = xx.*xx;
% 	NMDAgate = xx./(1+xx);
%     
%      
%     %g= ( Gampa(i+1)       +NMDAgate.*Gnmda(i+1)        + Ggabaa(i+1)         + Ggabab(i+1));
% 	%E= ( Gampa(i+1)*E_AMPA+NMDAgate.*Gnmda(i+1)*E_NMDA + Ggabaa(i+1)*E_GABAa + Ggabab(i+1)*E_GABAb);
%     %I(i) = v(i).*g-E;
%     if g_true_I_false
%         I(i) = (v(i)-E_AMPA)*Gampa(i+1) + (v(i)-E_NMDA)*NMDAgate*Gnmda(i+1)+ ...
%             (v(i)-E_GABAa)*Ggabaa(i+1)+(v(i)-E_GABAb)*Ggabab(i+1);
%     else
%         I(i) = se(i);
%     end
%     
%     v1=v(i); u1=u(i);
%     
%     v1=v1+0.5*(k.*(v1-vr).*(v1-vt)-u1-I(i))./C;    % for numerical 
%     %v(v>104)=104;   % Nothing higher than Eca2+
%     v1(v1>51)=51;   % Nothing higher than typical action potential
%     v1(v1<-90)=-90;  % cheating on numerical integration. Nothing lower than EK
% 
%     u1=u1+0.5*a.*(b.*(v1-vr)-u1);                   % step is 0.5 ms
%  
%     v1=v1+0.5*(k.*(v1-vr).*(v1-vt)-u1-I(i))./C;    % for numerical 
%     %v(v>104)=104;   % Nothing higher than Eca2+
%     v1(v1>51)=51;   % Nothing higher than Eca2+
%     v1(v1<-90)=-90;  % cheating on numerical integration. Nothing lower than EK
% 
%     u1=u1+0.5*a.*(b.*(v1-vr)-u1);                   % step is 0.5 ms
% 
%     v(i+1)=v1; u(i+1)=u1;
%     
%     %v(i+1)=v(i)+.5*dt_ms*(k*(v(i)-vr)*(v(i)-vt)-u(i)-I(i))/C;
%     %v(i+1)=v(i+1)+.5*dt_ms*(k*(v(i+1)-vr)*(v(i+1)-vt)-u(i)-I(i))/C;
%     %u(i+1)=u(i)+dt_ms*a*(b*(v(i)-vr)-u(i));
%     if v(i+1)>=vp
%         v(i)=vp;
%         v(i+1)=c;
%         u(i+1)=u(i+1)+d;
%         spikes = spikes + 1;
%     end
% end
% figure(1)
% subplot(3,1,1)
% plot(dt_ms*(1:n),[v' u']);
% subplot(3,1,2)
% plot(dt_ms*(1:n),[Gampa' Gnmda' Ggabaa' Ggabab']);
% title('conductances');
% subplot(3,1,3);
% plot(dt_ms*(1:n),I);
% title('Isyn (pA)');
% drawnow
